# --- Core Libraries ---
import pandas as pd
import numpy as np
import warnings
import os
import joblib
from typing import List, Dict, Tuple, Any

# --- Machine Learning ---
from sklearn.model_selection import TimeSeriesSplit
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD
from sklearn.pipeline import Pipeline
from sklearn.linear_model import MultiTaskLassoCV
import skl2onnx
from skl2onnx.common.data_types import FloatTensorType
import onnxruntime as rt

# --- Visualization & Utilities ---
import plotly.io as pio
from rich.console import Console
from rich.table import Table

# --- Custom Project Modules ---
import utils.custom_plotly_style_te as cpte
import utils.eda_visualizations as eda
import utils.analysis_helpers as ah
import utils.evaluations as eval
import utils.diagnostics as diag
import utils.training as train
import utils.drift_detection as dd
import utils.hourly_processing as hp

class PipelineConfig:
    """
    A centralized configuration class for managing paths, parameters, and flags
    throughout the notebook. Reproducibility? Yea. Quick and easy key settings
    modifications? Also yea:)
    """

    # --- Path & File Configuration ---
    DAILY_DATA_PATH: str = "data/weather_hcm_daily.csv"
    HOURLY_DATA_PATH = "data/weather_hcm_hourly.csv"
    UTILS_DIR: str = "utils"

    # --- General Model Parameters ---
    GLOBAL_RANDOM_SEED: int = 105
    HORIZONS: List[int] = [1, 2, 3, 4, 5] # Predict from t+1 to t+5
    TARGET_VARIABLE: str = "temp"

    # --- Data Split Parameters ---
    # Using an 85/15 split for a substantial test set. We use TimeSeriesSplit CV so can do without three-way holdout
    TEST_SIZE: float = 0.15


config = PipelineConfig()
os.makedirs(config.UTILS_DIR, exist_ok=True)
console = Console()

warnings.filterwarnings("ignore")

pio.templates["te"] = cpte.te_style_template
pio.templates.default = "te"
pd.options.plotting.backend = "plotly"
pio.renderers.default = "notebook+svg"
# svg_renderer = pio.renderers["svg"]
# svg_renderer.width = 1400
# svg_renderer.height = 900

# Load the dataset from the configured path
df_raw = pd.read_csv(config.DAILY_DATA_PATH)

# Convert 'datetime' column to datetime objects and set as the index
df_raw["datetime"] = pd.to_datetime(df_raw["datetime"])
df_raw.set_index("datetime", inplace=True)

# --- Perform a high-level data audit ---
console.print("--- [bold]DataFrame Info[/bold] ---")
df_raw.info(verbose=False, memory_usage="deep")

console.print("\n--- [bold]Descriptive Statistics (Numerical Features)[/bold] ---")
display(df_raw.describe().T)

console.print("\n--- [bold]Missing Value Percentages[/bold] ---")
missing_values = df_raw.isnull().sum()
missing_percentage = (missing_values / len(df_raw) * 100).round(2)
missing_df = pd.DataFrame({
    "Missing Count": missing_values,
    "Missing Percentage": missing_percentage,
}).sort_values(by="Missing Percentage", ascending=False)

console.print(missing_df[missing_df["Missing Count"] > 0])

--- DataFrame Info ---

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 3927 entries, 2015-01-01 to 2025-10-01
Columns: 37 entries, name to source
dtypes: float64(25), int64(2), object(10)
memory usage: 3.7 MB

--- Descriptive Statistics (Numerical Features) ---

--- Missing Value Percentages ---

            Missing Count  Missing Percentage
severerisk           2731               69.54
preciptype            924               23.53

fig_target_timeseries = eda.plot_target_variable_time_series(df=df_raw, target_col=config.TARGET_VARIABLE)
fig_target_timeseries.show()

fig_monthly_seasonality = eda.plot_monthly_seasonality(df=df_raw, target_col=config.TARGET_VARIABLE)
fig_monthly_seasonality.show()

# fig_stl = eda.plot_stl_decomposition(df=df_raw, target_col=config.TARGET_VARIABLE)
fig_stl = pio.read_json("./assets/plots/fig_stl_decomposition.json")
fig_stl.show()

# Select a representative subset of key numerical features - a mix of atmospheric, precipitation, and solar variables
features_to_plot = [
    "humidity",
    "dew",
    "sealevelpressure",
    "cloudcover",
    "solarradiation",
    "uvindex",
    "precip",
    "windspeed",
    "windgust",
]

feature_color_map = {
    feature: pio.templates["te"].layout.colorway[i % len(pio.templates["te"].layout.colorway)]
    for i, feature in enumerate(features_to_plot)
}

fig_distributions = eda.plot_feature_distributions(df=df_raw, features=features_to_plot, color_map=feature_color_map)
fig_distributions.show()

pairplot_features = [
    "temp",
    "dew",
    "precip",
    "humidity",
    "windspeed",
    "solarradiation",
    "cloudcover",
]

fig_pairplot = eda.plot_seasonal_pairplot(df=df_raw, features=pairplot_features)
fig_pairplot.show()

df_processed = df_raw.copy()

cols_to_drop = [
    # Constant Identifiers
    "name",
    "address",
    "resolvedAddress",
    "latitude",
    "longitude",
    "source",
    # Duplicate
    "windspeedmax",
    # Descriptive/Object columns to be engineered later
    "description",
    "conditions",
    "icon",
    "sunrise",
    "sunset",
    # Redundant energy column
    "solarenergy",
]

shape_before = df_processed.shape

df_processed.drop(columns=cols_to_drop, inplace=True, errors="ignore")

shape_after = df_processed.shape

console.print(f"DataFrame shape before sanitization: [yellow]{shape_before}[/yellow]")
console.print(f"DataFrame shape after sanitization:  [green]{shape_after}[/green]")
console.print(f"Columns removed: [red]{shape_before[1] - shape_after[1]}[/red]")

console.print("\n[bold]Remaining columns in the processed DataFrame:[/bold]")
console.print(df_processed.columns.tolist())

DataFrame shape before sanitization: (3927, 37)

DataFrame shape after sanitization:  (3927, 24)

Columns removed: 13

Remaining columns in the processed DataFrame:

[
    'tempmax',
    'tempmin',
    'temp',
    'feelslikemax',
    'feelslikemin',
    'feelslike',
    'dew',
    'humidity',
    'precip',
    'precipprob',
    'precipcover',
    'preciptype',
    'windgust',
    'windspeed',
    'windspeedmean',
    'windspeedmin',
    'winddir',
    'sealevelpressure',
    'cloudcover',
    'visibility',
    'solarradiation',
    'uvindex',
    'severerisk',
    'moonphase'
]

# 1. Fill 'preciptype' based on the deterministic relationship with 'precip'
df_processed["preciptype"].fillna("none", inplace=True)

# 2. Fill 'severerisk' with 0, as missingness implies no recorded severe risk
df_processed["severerisk"].fillna(0, inplace=True)


remaining_nulls = df_processed.isnull().sum().sum()

null_count_color = "green" if remaining_nulls == 0 else "red"
console.print(
    f"Total remaining null values in the DataFrame: [bold {null_count_color}]{remaining_nulls}[/bold {null_count_color}]"
)

if remaining_nulls == 0:
    console.print("\n[bold]Value counts for 'preciptype' after imputation:[/bold]")
    print(df_processed["preciptype"].value_counts())

    console.print("\n[bold]Value counts for 'severerisk' after imputation:[/bold]")
    print(df_processed["severerisk"].value_counts())

else:
    console.print("There are still missing values that need to be addressed.")

Total remaining null values in the DataFrame: 0

Value counts for 'preciptype' after imputation:

preciptype
rain    3003
none     924
Name: count, dtype: int64

Value counts for 'severerisk' after imputation:

severerisk
0.0     2731
30.0     546
10.0     448
60.0     195
75.0       6
5.0        1
Name: count, dtype: int64

df_numerical = df_processed.select_dtypes(include=np.number)

fig_pearson = eda.plot_correlation_heatmap_enhanced(df=df_numerical, method="pearson")
fig_pearson.show()

ah.display_correlation_summary(df=df_numerical, target=config.TARGET_VARIABLE, method="pearson")

Target Correlation Summary (Pearson) 

  Most Positively Correlated                            Most Negatively Correlated with                            
         with 'temp'                                                 'temp'                                        
┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓                         ┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓                          
┃ Feature      ┃ Coefficient ┃                         ┃ Feature          ┃ Coefficient ┃                          
┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩                         ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩                          
│ feelslike    │       0.893 │                         │ humidity         │      -0.340 │                          
│ tempmax      │       0.847 │                         │ cloudcover       │      -0.308 │                          
│ tempmin      │       0.792 │                         │ precip           │      -0.222 │                          
│ feelslikemin │       0.771 │                         │ sealevelpressure │      -0.128 │                          
│ feelslikemax │       0.634 │                         │ winddir          │      -0.087 │                          
└──────────────┴─────────────┘                         └──────────────────┴─────────────┘

Predictor Multicollinearity Summary (Pearson) 

       Highest Positive Predictor Pairs                           Highest Negative Predictor Pairs                 
┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓            ┏━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓         
┃ Feature 1      ┃ Feature 2    ┃ Coefficient ┃            ┃ Feature 1  ┃ Feature 2        ┃ Coefficient ┃         
┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩            ┡━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩         
│ tempmin        │ feelslikemin │       0.911 │            │ dew        │ sealevelpressure │      -0.652 │         
│ solarradiation │ uvindex      │       0.911 │            │ cloudcover │ solarradiation   │      -0.601 │         
│ temp           │ feelslike    │       0.893 │            │ humidity   │ solarradiation   │      -0.567 │         
│ tempmax        │ temp         │       0.847 │            │ cloudcover │ uvindex          │      -0.557 │         
│ feelslikemax   │ feelslike    │       0.833 │            │ humidity   │ sealevelpressure │      -0.554 │         
└────────────────┴──────────────┴─────────────┘            └────────────┴──────────────────┴─────────────┘

fig_spearman = eda.plot_correlation_heatmap_enhanced(df=df_numerical, method="spearman")
fig_spearman.show()

ah.display_correlation_summary(df=df_numerical, target=config.TARGET_VARIABLE, method="spearman")

Target Correlation Summary (Spearman) 

Most Positively Correlated with                           Most Negatively Correlated with                          
             'temp'                                                    'temp'                                      
┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓                         ┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓                        
┃ Feature        ┃ Coefficient ┃                         ┃ Feature          ┃ Coefficient ┃                        
┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩                         ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩                        
│ feelslike      │       0.874 │                         │ humidity         │      -0.393 │                        
│ tempmax        │       0.840 │                         │ cloudcover       │      -0.335 │                        
│ tempmin        │       0.757 │                         │ precip           │      -0.160 │                        
│ feelslikemin   │       0.756 │                         │ winddir          │      -0.101 │                        
│ solarradiation │       0.599 │                         │ sealevelpressure │      -0.091 │                        
└────────────────┴─────────────┘                         └──────────────────┴─────────────┘

Predictor Multicollinearity Summary (Spearman) 

       Highest Positive Predictor Pairs                           Highest Negative Predictor Pairs                 
┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓            ┏━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┓         
┃ Feature 1      ┃ Feature 2    ┃ Coefficient ┃            ┃ Feature 1  ┃ Feature 2        ┃ Coefficient ┃         
┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩            ┡━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━┩         
│ tempmin        │ feelslikemin │       0.999 │            │ dew        │ sealevelpressure │      -0.610 │         
│ solarradiation │ uvindex      │       0.881 │            │ humidity   │ solarradiation   │      -0.600 │         
│ temp           │ feelslike    │       0.874 │            │ cloudcover │ solarradiation   │      -0.584 │         
│ tempmax        │ temp         │       0.840 │            │ humidity   │ visibility       │      -0.566 │         
│ feelslikemax   │ feelslike    │       0.823 │            │ humidity   │ uvindex          │      -0.559 │         
└────────────────┴──────────────┴─────────────┘            └────────────┴──────────────────┴─────────────┘

console.print(
    f"[bold]Final Shape of Processed DataFrame:[/bold] [green]{df_processed.shape}[/green]",
    "\n[bold]Data Types and Memory Usage:[/bold]",
)
df_processed.info(memory_usage="deep")

Final Shape of Processed DataFrame: (3927, 24) 
Data Types and Memory Usage:

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 3927 entries, 2015-01-01 to 2025-10-01
Data columns (total 24 columns):
 #   Column            Non-Null Count  Dtype  
---  ------            --------------  -----  
 0   tempmax           3927 non-null   float64
 1   tempmin           3927 non-null   float64
 2   temp              3927 non-null   float64
 3   feelslikemax      3927 non-null   float64
 4   feelslikemin      3927 non-null   float64
 5   feelslike         3927 non-null   float64
 6   dew               3927 non-null   float64
 7   humidity          3927 non-null   float64
 8   precip            3927 non-null   float64
 9   precipprob        3927 non-null   int64  
 10  precipcover       3927 non-null   float64
 11  preciptype        3927 non-null   object 
 12  windgust          3927 non-null   float64
 13  windspeed         3927 non-null   float64
 14  windspeedmean     3927 non-null   float64
 15  windspeedmin      3927 non-null   float64
 16  winddir           3927 non-null   float64
 17  sealevelpressure  3927 non-null   float64
 18  cloudcover        3927 non-null   float64
 19  visibility        3927 non-null   float64
 20  solarradiation    3927 non-null   float64
 21  uvindex           3927 non-null   int64  
 22  severerisk        3927 non-null   float64
 23  moonphase         3927 non-null   float64
dtypes: float64(21), int64(2), object(1)
memory usage: 939.6 KB

df_engineered = df_processed.copy()

# --- 4.1.2: Basic Temporal Features ---
# These features help the model understand the long-term trend and basic monthly groupings.
df_engineered["year"] = df_engineered.index.year
df_engineered["month"] = df_engineered.index.month
df_engineered["week_of_year"] = df_engineered.index.isocalendar().week.astype(int)
df_engineered["day_of_year"] = df_engineered.index.dayofyear

# --- 4.1.3: Cyclical Features ---
# Sine/cosine transformations are essential for cyclical data. They encode the
# "closeness" of time points (e.g., Dec 31 is close to Jan 1) in a way that
# linear models can understand.
day_of_year_max = 366.0  # Account for leap years
month_max = 12.0

df_engineered["sin_day_of_year"] = np.sin(2 * np.pi * df_engineered["day_of_year"] / day_of_year_max)
df_engineered["cos_day_of_year"] = np.cos(2 * np.pi * df_engineered["day_of_year"] / day_of_year_max)
df_engineered["sin_month"] = np.sin(2 * np.pi * df_engineered["month"] / month_max)
df_engineered["cos_month"] = np.cos(2 * np.pi * df_engineered["month"] / month_max)

# --- 4.1.4: Fourier Series Features ---
# HCMC has a bimodal seasonal pattern. While the annual sin/cos
# captures the main cycle, higher-frequency Fourier terms can help model this
# more complex, multi-frequency seasonality (e.g., a semi-annual cycle). Just maybe, but hey it doesn't hurt to try.
for k in range(2, 4):  # Add k=2 (semi-annual) and k=3 (quarterly) terms
    df_engineered[f"fourier_sin_{k}"] = np.sin(2 * np.pi * k * df_engineered["day_of_year"] / day_of_year_max)
    df_engineered[f"fourier_cos_{k}"] = np.cos(2 * np.pi * k * df_engineered["day_of_year"] / day_of_year_max)


temporal_features_added = [
    "year",
    "month",
    "week_of_year",
    "day_of_year",
    "sin_day_of_year",
    "cos_day_of_year",
    "sin_month",
    "cos_month",
    "fourier_sin_2",
    "fourier_cos_2",
    "fourier_sin_3",
    "fourier_cos_3",
]
console.print(
    "[bold]Temporal features created.[/bold]\n",
    f" - Number of new features: [green]{len(temporal_features_added)}[/green]\n",
    f" - DataFrame shape is now: [green]{df_engineered.shape}[/green]",
)

console.print("\n[bold]Sample of new temporal features:[/bold]")
display(df_engineered[temporal_features_added].head())

Temporal features created.
  - Number of new features: 12
  - DataFrame shape is now: (3927, 36)

Sample of new temporal features:

# --- 4.2.1: Daylight Hours ---
# This is a more direct measure of solar energy potential than day_of_year.
sunrise_sunset = pd.read_csv(
    config.DAILY_DATA_PATH, usecols=["datetime", "sunrise", "sunset"], parse_dates=["datetime", "sunrise", "sunset"]
).set_index("datetime")

df_engineered["daylight_hours"] = (sunrise_sunset["sunset"] - sunrise_sunset["sunrise"]).dt.total_seconds() / 3600

# --- 4.2.2: Daily Temperature Range ---
# This feature quantifies daily temperature volatility. A low range often indicates
# a cloudy or rainy day, while a high range suggests a clear, sunny day.
df_engineered["temp_range"] = df_engineered["tempmax"] - df_engineered["tempmin"]

# --- 4.2.3: Dew Point Deficit ---
# This measures how close the air is to saturation. A smaller deficit means
# higher relative humidity and can inhibit nighttime cooling.
df_engineered["dew_point_deficit"] = df_engineered["temp"] - df_engineered["dew"]


domain_features_added = ["daylight_hours", "temp_range", "dew_point_deficit"]
console.print(
    "[bold]Domain-specific features created.[/bold]\n",
    f" - Number of new features: [green]{len(domain_features_added)}[/green]\n",
    f" - DataFrame shape is now: [green]{df_engineered.shape}[/green]",
)

console.print("\n[bold]Sample of new domain-specific features:[/bold]")
display(df_engineered[domain_features_added].head())

Domain-specific features created.
  - Number of new features: 3
  - DataFrame shape is now: (3927, 39)

Sample of new domain-specific features:

# --- 4.3.1: Lag Features ---
# Define the features to lag and the time periods for the lookback.
# We choose key dynamic variables that characterize the daily weather system.
features_to_lag = [
    "temp",
    "humidity",
    "solarradiation",
    "precip",
    "windspeed",
    "cloudcover",
    "dew",
    "sealevelpressure",
    "windgust",
]
lag_periods = [1, 2, 3, 7, 14, 28]

for feature in features_to_lag:
    for lag in lag_periods:
        df_engineered[f"{feature}_lag_{lag}"] = df_engineered[feature].shift(lag)

# --- 4.3.2: Rolling Window Features ---
# Define features and windows for rolling statistics.
features_for_rolling = [
    "temp",
    "humidity",
    "windspeed",
    "solarradiation",
    "dew",
    "cloudcover",  # Precip has its own logic
]
window_sizes = [7, 14, 28]
rolling_stats = ["mean", "std"]

for feature in features_for_rolling:
    # Shift the series by 1 to prevent data leakage.
    series_shifted = df_engineered[feature].shift(1)
    for window in window_sizes:
        for stat in rolling_stats:
            col_name = f"{feature}_roll_{window}d_{stat}"
            df_engineered[col_name] = series_shifted.rolling(window=window, min_periods=1).agg(stat)

# Specific logic for precipitation: sum is more meaningful than mean.
precip_shifted = df_engineered["precip"].shift(1)
for window in window_sizes:
    df_engineered[f"precip_roll_{window}d_sum"] = precip_shifted.rolling(window=window, min_periods=1).sum()


shape_before_fe = df_processed.shape[1]
shape_after_fe = df_engineered.shape[1]
num_new_features = shape_after_fe - shape_before_fe - len(temporal_features_added) - len(domain_features_added)

console.print(
    "[bold]Lag and rolling features created.[/bold]\n",
    f" - Number of new lag/rolling features: [green]{num_new_features}[/green]\n",
    f" - DataFrame shape is now: [green]{df_engineered.shape}[/green]",
)

Lag and rolling features created.
  - Number of new lag/rolling features: 93
  - DataFrame shape is now: (3927, 132)

# Perform One-Hot Encoding on the 'preciptype' column
df_engineered = pd.get_dummies(
    df_engineered,
    columns=["preciptype"],
    prefix="preciptype",
    dtype=int,
)


encoded_cols = [col for col in df_engineered.columns if "preciptype_" in col]
console.print(
    "[bold]Original 'preciptype' column removed.[/bold]\n",
    f" - New binary columns added: `{encoded_cols}`\n",
    f" - DataFrame shape is now: [green]{df_engineered.shape}[/green]",
)

console.print("\n[bold]Sample of new encoded features:[/bold]")
display(df_engineered[encoded_cols].head())

Original 'preciptype' column removed.
  - New binary columns added: `['preciptype_none', 'preciptype_rain']`
  - DataFrame shape is now: (3927, 133)

Sample of new encoded features:

text_data = df_raw[["description"]].copy()
df_engineered = df_engineered.join(text_data)

# --- Define the NLP processing pipeline ---
n_svd_components = 7  # The `description` field has 26 distinct features, mostly repetitive, so 7 is sufficient, I think
nlp_pipeline = Pipeline([
    (
        "tfidf",
        TfidfVectorizer(
            stop_words="english",
            ngram_range=(1, 2),  # Capture both single words and two-word phrases
            max_features=1000,  # Limit vocabulary to the top 1000 terms
        ),
    ),
    ("svd", TruncatedSVD(n_components=n_svd_components, random_state=config.GLOBAL_RANDOM_SEED)),
])

text_features = nlp_pipeline.fit_transform(df_engineered["description"])

text_feature_names = [f"text_svd_{i}" for i in range(n_svd_components)]
df_text_features = pd.DataFrame(text_features, index=df_engineered.index, columns=text_feature_names)

df_engineered = pd.concat([df_engineered, df_text_features], axis=1)

df_engineered.drop(columns=["description"], inplace=True)


console.print(
    f" - Number of new text-based features (SVD components): [green]{n_svd_components}[/green]\n",
    f" - DataFrame shape is now: [green]{df_engineered.shape}[/green]",
    sep="",
)

console.print("\n[bold]Sample of new text features:[/bold]")
display(df_engineered[text_feature_names].head())

 - Number of new text-based features (SVD components): 7
 - DataFrame shape is now: (3927, 140)

Sample of new text features:

max_lag = max(lag_periods)
max_rolling_window = max(window_sizes)
max_lookback = max(max_lag, max_rolling_window)

console.print(
    f"Longest lag period: [cyan]{max_lag}[/cyan] days\n",
    f"Longest rolling window: [cyan]{max_rolling_window}[/cyan] days\n",
    f"==> Maximum lookback period: [bold yellow]{max_lookback}[/bold yellow] days\n",
    f"\nThis means the first {max_lookback} rows will contain NaN values and will be dropped during the modeling prep stage.\n",
    f"[bold]Final Shape of DataFrame:[/bold] [green]{df_engineered.shape}[/green]",
    sep="",
)

console.print("[bold]First 5 Rows (showing initial NaNs):[/bold]")
display(df_engineered.head())

Longest lag period: 28 days
Longest rolling window: 28 days
==> Maximum lookback period: 28 days

This means the first 28 rows will contain NaN values and will be dropped during the modeling prep stage.
Final Shape of DataFrame: (3927, 140)

First 5 Rows (showing initial NaNs):

df_model_ready = df_engineered.copy()

# --- 1. Create the Multi-Horizon Target Variables ---
for h in config.HORIZONS:
    df_model_ready[f"target_temp_t+{h}"] = df_model_ready["temp"].shift(-h)

# --- 2. Final Alignment via NaN Removal ---
shape_before_align = df_model_ready.shape
df_model_ready.dropna(inplace=True)
shape_after_align = df_model_ready.shape

rows_dropped = shape_before_align[0] - shape_after_align[0]

console.print(
    f"Shape before final alignment: [green]{shape_before_align}[/green]\n",
    f"Shape after final alignment:  [green]{shape_after_align}[/green]\n",
    f"Rows dropped due to lookback/lookahead: [red]{rows_dropped}[/red]",
    sep="",
)

Shape before final alignment: (3927, 145)
Shape after final alignment:  (3894, 145)
Rows dropped due to lookback/lookahead: 33

# --- 1. Define Feature and Target Columns ---
target_cols = [f"target_temp_t+{h}" for h in config.HORIZONS]

# We must drop the original 'temp' column as is now a source of data leakage (because it's
# from the same day as the target values we are trying to predict)
features_to_drop = target_cols + ["temp"]

# All remaining columns are predictors.
feature_cols = [col for col in df_model_ready.columns if col not in features_to_drop]

X = df_model_ready[feature_cols]
y = df_model_ready[target_cols]

# --- 2. Perform the Chronological Split ---
n_total = len(X)
n_test = int(n_total * config.TEST_SIZE)
n_train = n_total - n_test

X_train_full = X.iloc[:n_train]
y_train = y.iloc[:n_train]

X_test_full = X.iloc[n_train:]
y_test = y.iloc[n_train:]


console.print(
    f" - Total samples:     {n_total}\n",
    f" - Training set size: {n_train} ({1 - config.TEST_SIZE:.0%})\n",
    f" - Test set size:     {n_test} ({config.TEST_SIZE:.0%})",
    sep="",
)

console.rule("[bold]Time Range of the sets:[/bold]", align="center")
console.print(
    f"Training Set Time Range:  {X_train_full.index.min().date()} to {X_train_full.index.max().date()}\n",
    f"Test Set Time Range:      {X_test_full.index.min().date()} to {X_test_full.index.max().date()}",
    sep="",
)

 - Total samples:     3894
 - Training set size: 3310 (85%)
 - Test set size:     584 (15%)

───────────────────────────────────────────── Time Range of the sets: ─────────────────────────────────────────────

Training Set Time Range:  2015-01-29 to 2024-02-20
Test Set Time Range:      2024-02-21 to 2025-09-26

X_train_tree = X_train_full
X_test_tree = X_test_full

# --- Use MultiTaskLassoCV to select a low-collinearity subset for linear models ---
# Fit only on the training data
# lasso_pipeline = Pipeline([
#     ("scaler", RobustScaler()),
#     (
#         "lasso",
#         MultiTaskLassoCV(
#             cv=TimeSeriesSplit(n_splits=5),
#             random_state=config.GLOBAL_RANDOM_SEED,
#             n_jobs=-1,
#         ),
#     ),
# ])

# lasso_pipeline.fit(X_train_full, y_train)

# # Identify features with non-zero coefficients.
# coefficients = lasso_pipeline.named_steps["lasso"].coef_
# selected_features_mask = np.abs(coefficients).sum(axis=0) > 0
selected_features_mask = joblib.load("./assets/saved_results/selected_features_mask.pkl") # Load cho nhanh
linear_feature_cols = X_train_full.columns[selected_features_mask].tolist()

X_train_linear = X_train_full[linear_feature_cols]
X_test_linear = X_test_full[linear_feature_cols]


console.print(
    f" - Tree-based model feature count: [cyan]{X_train_tree.shape[1]}[/cyan]\n",
    f" - Linear model feature count:     [cyan]{X_train_linear.shape[1]}[/cyan] (selected by LASSO)",
    sep="",
)

 - Tree-based model feature count: 139
 - Linear model feature count:     41 (selected by LASSO)

ah.display_lasso_summary_table(
    all_features=X_train_tree.columns.tolist(),
    selected_features=X_train_linear.columns.tolist(),
)

fig_sunburst_final = ah.plot_lasso_summary_sunburst(
    all_features=X_train_tree.columns.tolist(),
    selected_features=X_train_linear.columns.tolist(),
)
fig_sunburst_final.show()

           LASSO Feature Selection Summary           
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━┳━━━━━━━━━┳━━━━━━━┓
┃ Feature Category         ┃ Kept ┃ Dropped ┃ Total ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━╇━━━━━━━━━╇━━━━━━━┩
│ Raw Predictor            │  13  │    9    │  22   │
│ Temporal (Cyclical)      │  4   │    4    │   8   │
│ Text (SVD)               │  3   │    4    │   7   │
│ Lag (Humidity)           │  1   │    5    │   6   │
│ Lag (Dew)                │  2   │    4    │   6   │
│ Lag (Cloudcover)         │  1   │    5    │   6   │
│ Rolling (Dew)            │  0   │    6    │   6   │
│ Rolling (Cloudcover)     │  2   │    4    │   6   │
│ Lag (Precip)             │  6   │    0    │   6   │
│ Lag (Sealevelpressure)   │  2   │    4    │   6   │
│ Lag (Windgust)           │  0   │    6    │   6   │
│ Lag (Windspeed)          │  0   │    6    │   6   │
│ Lag (Temp)               │  3   │    3    │   6   │
│ Lag (Solarradiation)     │  0   │    6    │   6   │
│ Rolling (Windspeed)      │  1   │    5    │   6   │
│ Rolling (Temp)           │  2   │    4    │   6   │
│ Rolling (Humidity)       │  0   │    6    │   6   │
│ Rolling (Solarradiation) │  1   │    5    │   6   │
│ Temporal (Basic)         │  0   │    4    │   4   │
│ Domain-Specific          │  0   │    3    │   3   │
│ Rolling (Precip)         │  0   │    3    │   3   │
│ Categorical              │  0   │    2    │   2   │
└──────────────────────────┴──────┴─────────┴───────┘

tscv_splitter = TimeSeriesSplit(n_splits=5)

baseline_cv_scores = eval.evaluate_persistence_model(
    df_model_ready=df_model_ready,
    X_train=X_train_full,
    y_train_all_horizons=y_train,
    horizons=config.HORIZONS,
    cv_splitter=tscv_splitter,
)

baseline_cv_scores_for_table = {
    horizon: data["scores"]["r2"] for horizon, data in baseline_cv_scores["Baseline (Persistence)"].items()
}
eval.display_baseline_results(baseline_cv_scores_for_table)

   Persistence (Naive) Baseline Performance    
       (Cross-Validated on Training Set)       
┏━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┓
┃ Horizon ┃ Mean CV R² ┃ Std Dev of CV Scores ┃
┡━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━┩
│   t+1   │   0.5455   │        0.0675        │
│   t+2   │   0.2080   │        0.1193        │
│   t+3   │   0.0634   │        0.1463        │
│   t+4   │   0.0172   │        0.1399        │
│   t+5   │  -0.0467   │        0.1391        │
│ Average │   0.1575   │        0.1224        │
└─────────┴────────────┴──────────────────────┘

results_files_to_load = [
    "./assets/saved_results/bakeoff_results_linear.pkl",
    "./assets/saved_results/bakeoff_results_tree.pkl",
    "./assets/saved_results/bakeoff_results_advanced_ts.pkl",
    "./assets/saved_results/bakeoff_results_dl.pkl",
]

df_leaderboard = eval.generate_leaderboard_data(results_files=results_files_to_load, baseline_scores=baseline_cv_scores)

eval.display_leaderboard(df_leaderboard, title="Stage 1 - High-level Leaderboard")

finalist_models = [
    "RidgeCV",
    "LinearRegression",
    "ElasticNetCV",
    "CatBoostRegressor",
    "LGBMRegressor",
]

                                         Stage 1 - High-level Leaderboard                                          
                               (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┓
┃      ┃                         ┃                 ┃             ┃               ┃              ┃ Avg CV MAPE (%) ┃
┃ Rank ┃ Model                   ┃  Model Family   ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃        ↓        ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━┩
│  #1  │ RidgeCV                 │     Linear      │   0.4219    │    0.9570     │    0.7514    │      2.68%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #2  │ LinearRegression        │     Linear      │   0.4185    │    0.9595     │    0.7523    │      2.68%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #3  │ ElasticNetCV            │     Linear      │   0.4088    │    0.9673     │    0.7596    │      2.71%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #4  │ CatBoostRegressor       │   Tree-based    │   0.3248    │    1.0258     │    0.8122    │      2.89%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #5  │ LGBMRegressor           │   Tree-based    │   0.2901    │    1.0485     │    0.8319    │      2.97%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #6  │ XGBRegressor            │   Tree-based    │   0.2370    │    1.0894     │    0.8629    │      3.07%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #7  │ Baseline (Persistence)  │    Baseline     │   0.1575    │    1.1479     │    0.8834    │      3.15%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #8  │ Prophet                 │  Advanced time  │   -1.2556   │    1.6671     │    1.3497    │      4.82%      │
│      │                         │     series      │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #9  │ LSTM                    │  Deep learning  │   -1.4927   │    1.9124     │    1.5530    │      5.46%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #10  │ SARIMAX                 │  Advanced time  │ -3118.7928  │    50.6686    │   44.6821    │     158.28%     │
│      │                         │     series      │             │               │              │                 │
└──────┴─────────────────────────┴─────────────────┴─────────────┴───────────────┴──────────────┴─────────────────┘

df_horizon_data = eval.generate_finalist_horizon_data(
    results_files=results_files_to_load,
    finalists=finalist_models,
    metrics=["r2", "rmse", "mae"],
)
eval.display_finalist_deep_dive_plots(df_horizon_data)

df_fold_data = eval.generate_finalist_fold_data(
    results_files=results_files_to_load,
    finalists=["RidgeCV", "CatBoostRegressor"],
    horizons_to_plot=["t+1", "t+5"],
)

fig_fold_plot = eval.display_finalist_fold_plot(df_fold_data)
fig_fold_plot.show()

results_files_to_load = [
    "./assets/saved_results/bakeoff_results_linear.pkl",
    "./assets/saved_results/bakeoff_results_tree.pkl",
    "./assets/saved_results/bakeoff_results_ensemble.pkl",
]

df_leaderboard_final = eval.generate_leaderboard_data(
    results_files=results_files_to_load,
    baseline_scores=baseline_cv_scores,
)
eval.display_leaderboard(df_leaderboard_final, title="Stage 3 - The Rise of Ensembles?")

                                         Stage 3 - The Rise of Ensembles?                                          
                               (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┓
┃      ┃                         ┃                ┃             ┃               ┃              ┃ Avg CV MAPE (%)  ┃
┃ Rank ┃ Model                   ┃  Model Family  ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃        ↓         ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━┩
│  #1  │ RidgeCV                 │     Linear     │   0.4219    │    0.9570     │    0.7514    │      2.68%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #2  │ LinearRegression        │     Linear     │   0.4185    │    0.9595     │    0.7523    │      2.68%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #3  │ SimpleAveragingEnsemble │ Ensemble (OOB) │   0.4175    │    0.9587     │    0.7540    │      2.69%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #4  │ ElasticNetCV            │     Linear     │   0.4088    │    0.9673     │    0.7596    │      2.71%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #5  │ ElasticNet_CatBoost_wit │ Ensemble (OOB) │   0.3839    │    0.9819     │    0.7741    │      2.76%       │
│      │ h_ElasticNetCV          │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #6  │ ElasticNetCV_LGBM_with_ │ Ensemble (OOB) │   0.3832    │    0.9819     │    0.7730    │      2.76%       │
│      │ RidgeCV                 │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #7  │ ElasticNet_CatBoost_wit │ Ensemble (OOB) │   0.3814    │    0.9836     │    0.7749    │      2.76%       │
│      │ h_RidgeCV               │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #8  │ RidgeCV_CatBoost_with_E │ Ensemble (OOB) │   0.3650    │    0.9931     │    0.7846    │      2.80%       │
│      │ lasticNeCV              │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #9  │ RidgeCV_CatBoost_with_R │ Ensemble (OOB) │   0.3627    │    0.9940     │    0.7853    │      2.80%       │
│      │ idgeCV                  │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│ #10  │ CatBoostRegressor       │   Tree-based   │   0.3248    │    1.0258     │    0.8122    │      2.89%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│ #11  │ LGBMRegressor           │   Tree-based   │   0.2901    │    1.0485     │    0.8319    │      2.97%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│ #12  │ XGBRegressor            │   Tree-based   │   0.2370    │    1.0894     │    0.8629    │      3.07%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│ #13  │ Baseline (Persistence)  │    Baseline    │   0.1575    │    1.1479     │    0.8834    │      3.15%       │
└──────┴─────────────────────────┴────────────────┴─────────────┴───────────────┴──────────────┴──────────────────┘

df_fold_data = eval.generate_finalist_fold_data(
    results_files=results_files_to_load,
    finalists=["RidgeCV", "SimpleAveragingEnsemble", "ElasticNet_CatBoost_with_ElasticNetCV"],
    horizons_to_plot=["t+1", "t+2", "t+3", "t+4", "t+5"],
)

fig_fold_plot = eval.display_finalist_fold_plot(df_fold_data)
fig_fold_plot.show()

results_files_to_load = [
    "./assets/saved_results/bakeoff_results_tree.pkl",
    "./assets/saved_results/evaluation_results_tuned_single_models_v2_weighted_CV.pkl",
]

df_leaderboard_final = eval.generate_leaderboard_data(results_files=results_files_to_load)
eval.display_leaderboard(df_leaderboard_final, title="Optuna on Single Models")

                                             Optuna on Single Models                                              
                              (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┓
┃ Rank ┃ Model                   ┃ Model Family ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃ Avg CV MAPE (%) ↓ ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━┩
│  #1  │ Tuned_LGBM_v2_weighted_ │  Tuned tree  │   0.4036    │    0.9674     │    0.7625    │       2.72%       │
│      │ CV                      │              │             │               │              │                   │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #2  │ Tuned_CatBoost_v2_weigh │  Tuned tree  │   0.3835    │    0.9822     │    0.7761    │       2.77%       │
│      │ ted_CV                  │              │             │               │              │                   │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #3  │ CatBoostRegressor       │  Tree-based  │   0.3248    │    1.0258     │    0.8122    │       2.89%       │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #4  │ LGBMRegressor           │  Tree-based  │   0.2901    │    1.0485     │    0.8319    │       2.97%       │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #5  │ XGBRegressor            │  Tree-based  │   0.2370    │    1.0894     │    0.8629    │       3.07%       │
└──────┴─────────────────────────┴──────────────┴─────────────┴───────────────┴──────────────┴───────────────────┘

results_files_to_load = [
    "./assets/saved_results/evaluation_results_tuned_single_models_v1_mean.pkl",
    "./assets/saved_results/evaluation_results_tuned_single_models_v2_weighted_CV.pkl",
    "./assets/saved_results/evaluation_results_tuned_LGBM_v3_linear_tree.pkl"
]

df_leaderboard_final = eval.generate_leaderboard_data(results_files=results_files_to_load)
eval.display_leaderboard(df_leaderboard_final, title="Average CV vs. Weighted CV optimization objective")

                                Average CV vs. Weighted CV optimization objective                                 
                              (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┓
┃ Rank ┃ Model                   ┃ Model Family ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃ Avg CV MAPE (%) ↓ ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━┩
│  #1  │ Tuned_LGBM_v2_weighted_ │  Tuned tree  │   0.4036    │    0.9674     │    0.7625    │       2.72%       │
│      │ CV                      │              │             │               │              │                   │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #2  │ Tuned_LGBM_v3_linear_tr │  Tuned tree  │   0.3938    │    0.9773     │    0.7674    │       2.74%       │
│      │ ee                      │              │             │               │              │                   │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #3  │ Tuned_LGBM_v1_mean      │  Tuned tree  │   0.3924    │    0.9749     │    0.7693    │       2.75%       │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #4  │ Tuned_CatBoost_v2_weigh │  Tuned tree  │   0.3835    │    0.9822     │    0.7761    │       2.77%       │
│      │ ted_CV                  │              │             │               │              │                   │
├──────┼─────────────────────────┼──────────────┼─────────────┼───────────────┼──────────────┼───────────────────┤
│  #5  │ Tuned_CatBoost_v1_mean  │  Tuned tree  │   0.3747    │    0.9901     │    0.7827    │       2.79%       │
└──────┴─────────────────────────┴──────────────┴─────────────┴───────────────┴──────────────┴───────────────────┘

results_files_to_load = [
    "./assets/saved_results/bakeoff_results_linear.pkl",
    "./assets/saved_results/bakeoff_results_tree.pkl",
    "./assets/saved_results/evaluation_results_tuned_single_models_v2_weighted_CV.pkl",
    "./assets/saved_results/bakeoff_results_ensemble.pkl",
]

df_fold_data = eval.generate_finalist_fold_data(
    results_files=results_files_to_load,
    finalists=["RidgeCV", "SimpleAveragingEnsemble", "CatBoostRegressor", "Tuned_LGBM_v2_weighted_CV"],
    horizons_to_plot=["t+1", "t+5"],
)

fig_fold_plot = eval.display_finalist_fold_plot(df_fold_data)
fig_fold_plot.show()

results_files_to_load = [
    "./assets/saved_results/bakeoff_results_linear.pkl",
    "./assets/saved_results/bakeoff_results_ensemble.pkl",
    "./assets/saved_results/evaluation_results_tuned_ensembles_test.pkl",
]

df_leaderboard_final = eval.generate_leaderboard_data(results_files=results_files_to_load)
eval.display_leaderboard(
    df_leaderboard_final,
    title="Ensembles - OOB vs Optuna",
    highlight_models=["Avg_Ridge_TunedLGBM", "SimpleAveragingEnsemble", "Stack_ElasticNet_LGBM_Ridge_Tuned", "RidgeCV"],
    exclude_models=["LinearRegression", "ElasticNetCV"],
)

                                             Ensembles - OOB vs Optuna                                             
                               (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┓
┃      ┃                         ┃                 ┃             ┃               ┃              ┃ Avg CV MAPE (%) ┃
┃ Rank ┃ Model                   ┃  Model Family   ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃        ↓        ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━┩
│  #1  │ Avg_Ridge_TunedLGBM     │  Tuned simple   │   0.4371    │    0.9432     │    0.7402    │      2.64%      │
│      │                         │       avg       │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #2  │ Avg_Ridge_TunedCatBoost │  Tuned simple   │   0.4331    │    0.9463     │    0.7437    │      2.65%      │
│      │                         │       avg       │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #3  │ Avg_ElasticNet_TunedLGB │  Tuned simple   │   0.4287    │    0.9498     │    0.7468    │      2.66%      │
│      │ M                       │       avg       │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #4  │ Avg_ElasticNet_TunedCat │  Tuned simple   │   0.4224    │    0.9546     │    0.7515    │      2.68%      │
│      │ Boost                   │       avg       │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #5  │ RidgeCV                 │     Linear      │   0.4219    │    0.9570     │    0.7514    │      2.68%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #7  │ SimpleAveragingEnsemble │ Ensemble (OOB)  │   0.4175    │    0.9587     │    0.7540    │      2.69%      │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #9  │ Stack_ElasticNet_LGBM_R │ Tuned stacking  │   0.4054    │    0.9663     │    0.7613    │      2.71%      │
│      │ idge_Tuned              │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #10  │ Stack_ElasticNet_CatBoo │ Tuned stacking  │   0.3926    │    0.9757     │    0.7688    │      2.74%      │
│      │ st_Ridge_Tuned          │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #11  │ ElasticNet_CatBoost_wit │ Ensemble (OOB)  │   0.3839    │    0.9819     │    0.7741    │      2.76%      │
│      │ h_ElasticNetCV          │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #12  │ ElasticNetCV_LGBM_with_ │ Ensemble (OOB)  │   0.3832    │    0.9819     │    0.7730    │      2.76%      │
│      │ RidgeCV                 │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #13  │ ElasticNet_CatBoost_wit │ Ensemble (OOB)  │   0.3814    │    0.9836     │    0.7749    │      2.76%      │
│      │ h_RidgeCV               │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #14  │ RidgeCV_CatBoost_with_E │ Ensemble (OOB)  │   0.3650    │    0.9931     │    0.7846    │      2.80%      │
│      │ lasticNeCV              │                 │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│ #15  │ RidgeCV_CatBoost_with_R │ Ensemble (OOB)  │   0.3627    │    0.9940     │    0.7853    │      2.80%      │
│      │ idgeCV                  │                 │             │               │              │                 │
└──────┴─────────────────────────┴─────────────────┴─────────────┴───────────────┴──────────────┴─────────────────┘

results_files_to_load = [
    "./assets/saved_results/evaluation_results_tuned_ensembles_test.pkl",
    "./assets/saved_results/evaluation_results_tuned_ensemble_weighted_avg.pkl",
]

df_leaderboard_final = eval.generate_leaderboard_data(results_files=results_files_to_load)
eval.display_leaderboard(
    df_leaderboard_final,
    title="Simple 50/50 vs Static Weighted-horizon Ensembles",
    exclude_models=[
        "Avg_Ridge_TunedCatBoost",
        "Avg_ElasticNet_TunedLGBM",
        "Avg_ElasticNet_TunedCatBoost",
        "Stack_ElasticNet_LGBM_Ridge_Tuned",
        "Stack_ElasticNet_CatBoost_Ridge_Tuned",
    ],
)

                                 Simple 50/50 vs Static Weighted-horizon Ensembles                                 
                               (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┓
┃      ┃                         ┃                 ┃             ┃               ┃              ┃ Avg CV MAPE (%) ┃
┃ Rank ┃ Model                   ┃  Model Family   ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃        ↓        ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━┩
│  #1  │ WeightedAvg_Ridge_Tuned │  Tuned simple   │   0.4416    │    0.9398     │    0.7405    │      2.64%      │
│      │ LGBM                    │       avg       │             │               │              │                 │
├──────┼─────────────────────────┼─────────────────┼─────────────┼───────────────┼──────────────┼─────────────────┤
│  #2  │ Avg_Ridge_TunedLGBM     │  Tuned simple   │   0.4371    │    0.9432     │    0.7402    │      2.64%      │
│      │                         │       avg       │             │               │              │                 │
└──────┴─────────────────────────┴─────────────────┴─────────────┴───────────────┴──────────────┴─────────────────┘

# df_fold_data = eval.generate_finalist_fold_data(
#     results_files=results_files_to_load,
#     finalists=["WeightedAvg_Ridge_TunedLGBM", "Avg_Ridge_TunedLGBM"],
#     horizons_to_plot=["t+1", "t+2", "t+3", "t+4", "t+5"],
# )

# fig_fold_plot = eval.display_finalist_fold_plot(df_fold_data)
# fig_fold_plot.show()

df_importance_lgbm = joblib.load("./assets/saved_results/df_importance_lgbm.pkl")

fig_bump_chart = eval.plot_importance_bump_chart(df_importance_lgbm)
fig_bump_chart.show()

diagnostic_data = joblib.load("./assets/saved_results/chosen_models_diagnostic_data.pkl")

# fig_resid_vs_time = diag.plot_residuals_vs_time(
#     diagnostic_data=diagnostic_data,
#     horizons_to_plot=["t+1", "t+5"],
# )
fig_resid_vs_time = pio.read_json("./assets/plots/fig_resid_vs_time.json")
fig_resid_vs_time.show()

fig_acf_pacf = diag.plot_acf_pacf(
    diagnostic_data=diagnostic_data,
    horizons_to_plot=["t+1"],
)

fig_acf_pacf.show()

fig_resid_vs_pred = diag.plot_residuals_vs_variable(
    diagnostic_data=diagnostic_data,
    horizon_to_plot="t+1",
    x_variable="predicted",
)
fig_resid_vs_pred.show()

fig_resid_vs_variable = diag.plot_residuals_vs_variable(
    diagnostic_data=diagnostic_data,
    horizon_to_plot="t+1",
    x_variable="solarradiation",
)
fig_resid_vs_variable.show()

shap_vs_perm_data = joblib.load("assets/saved_results/shap_vs_perm_data.pkl")

fig_shap_vs_perm = diag.plot_shap_vs_perm_comparison(
    shap_vs_perm_data=shap_vs_perm_data,
    horizons_to_plot=["t+1", "t+5"],
    n_top_features=20,
    sort_by="shap",
)
fig_shap_vs_perm.show()

df_enhancement_leaderboard = eval.generate_leaderboard_data(
    results_files=["./assets/saved_results/evaluation_results_ensemble_enhancements.pkl"]
)
eval.display_leaderboard(
    df_leaderboard=df_enhancement_leaderboard,
    title="Champion Ensemble: Base vs Enhanced",
)

                                        Champion Ensemble: Base vs Enhanced                                        
                               (Metrics averaged across all 5 horizons & 5 CV folds)                               
┏━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┓
┃      ┃                         ┃                ┃             ┃               ┃              ┃ Avg CV MAPE (%)  ┃
┃ Rank ┃ Model                   ┃  Model Family  ┃ Avg CV R² ↑ ┃ Avg CV RMSE ↓ ┃ Avg CV MAE ↓ ┃        ↓         ┃
┡━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━┩
│  #1  │ Champion (Unenhanced)   │    Unknown     │   0.4416    │    0.9398     │    0.7405    │      2.64%       │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #2  │ Ensemble + Isotonic     │ Ensemble (OOB) │   0.4365    │    0.9438     │    0.7406    │      2.64%       │
│      │ Calibration             │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #3  │ Ensemble + Residual     │ Tuned stacking │   0.4250    │    0.9524     │    0.7471    │      2.66%       │
│      │ Stacking (Tuned LGBM)   │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #4  │ Ensemble + Residual     │ Tuned stacking │   0.4249    │    0.9525     │    0.7472    │      2.66%       │
│      │ Stacking (Tuned LGBM) + │                │             │               │              │                  │
│      │ Calibration             │                │             │               │              │                  │
├──────┼─────────────────────────┼────────────────┼─────────────┼───────────────┼──────────────┼──────────────────┤
│  #5  │ Ensemble + Residual     │ Ensemble (OOB) │   0.1229    │    1.1290     │    0.8985    │      3.21%       │
│      │ Stacking (Ridge)        │                │             │               │              │                  │
└──────┴─────────────────────────┴────────────────┴─────────────┴───────────────┴──────────────┴──────────────────┘

final_models_to_evaluate = {
    "WeightedAvg_Ridge_TunedLGBM": {
        "type": "weighted_averaging_tuned",
    },
    "RidgeCV": {
        "type": "oob",
        "feature_set": "linear",
    },
}

final_predictions = train.load_and_predict_on_test(
    models_to_predict=final_models_to_evaluate,
    X_test_linear=X_test_linear,
    X_test_tree=X_test_tree,
    horizons=config.HORIZONS,
    load_dir="./assets/models/daily",
)

eval.display_final_performance_table(
    predictions=final_predictions,
    y_test=y_test,
    display_names={
        "WeightedAvg_Ridge_TunedLGBM": "Ensemble",
        "RidgeCV": "RidgeCV",
    },
)

                           Final Model Performance on Unseen Test Set                            
┏━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━┓
┃         ┃ Ensemble ┃ RidgeCV ┃ Ensemble ┃ RidgeCV  ┃ Ensemble ┃ RidgeCV ┃ Ensemble ┃ RidgeCV  ┃
┃ Horizon ┃  (R² ↑)  ┃ (R² ↑)  ┃ (RMSE ↓) ┃ (RMSE ↓) ┃ (MAE ↓)  ┃ (MAE ↓) ┃ (MAPE ↓) ┃ (MAPE ↓) ┃
┡━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━┩
│   t+1   │  0.7460  │ 0.7396  │  0.7912  │  0.8010  │  0.6265  │ 0.6324  │  2.19%   │  2.21%   │
├─────────┼──────────┼─────────┼──────────┼──────────┼──────────┼─────────┼──────────┼──────────┤
│   t+2   │  0.6306  │ 0.6232  │  0.9555  │  0.9651  │  0.7628  │ 0.7725  │  2.66%   │  2.71%   │
├─────────┼──────────┼─────────┼──────────┼──────────┼──────────┼─────────┼──────────┼──────────┤
│   t+3   │  0.5805  │ 0.5748  │  1.0195  │  1.0264  │  0.8130  │ 0.8144  │  2.84%   │  2.85%   │
├─────────┼──────────┼─────────┼──────────┼──────────┼──────────┼─────────┼──────────┼──────────┤
│   t+4   │  0.5689  │ 0.5658  │  1.0336  │  1.0374  │  0.8225  │ 0.8231  │  2.87%   │  2.88%   │
├─────────┼──────────┼─────────┼──────────┼──────────┼──────────┼─────────┼──────────┼──────────┤
│   t+5   │  0.5711  │ 0.5648  │  1.0325  │  1.0401  │  0.8246  │ 0.8299  │  2.88%   │  2.90%   │
├─────────┼──────────┼─────────┼──────────┼──────────┼──────────┼─────────┼──────────┼──────────┤
│ Average │  0.6194  │ 0.6136  │  0.9665  │  0.9740  │  0.7699  │ 0.7745  │  2.69%   │  2.71%   │
└─────────┴──────────┴─────────┴──────────┴──────────┴──────────┴─────────┴──────────┴──────────┘

model_name = "WeightedAvg_Ridge_TunedLGBM"
horizon_to_plot = "t+1"

y_true_plot = y_test[f"target_temp_{horizon_to_plot}"]
y_pred_plot = final_predictions[model_name][horizon_to_plot]

fig_ts_dashboard = diag.plot_performance_timeseries_dashboard(
    y_true=y_true_plot,
    y_pred=y_pred_plot,
    model_name="Champion Ensemble",
    horizon=horizon_to_plot,
    # add_rolling_mean=7,
)
fig_ts_dashboard.show()

fig_stat_dashboard = diag.plot_performance_statistical_dashboard(
    y_true=y_true_plot,
    y_pred=y_pred_plot,
    model_name="Champion Ensemble",
    horizon=horizon_to_plot,
)

fig_stat_dashboard.show()

# Load pre-trained ensemble models for each horizon
sim_ensemble = {f"h{i}": joblib.load(f"assets/models/daily/WeightedAvg_Ridge_TunedLGBM_h{i}.pkl") for i in range(1, 6)}

# Use the model-ready data that already has all engineered features
df_sim = df_model_ready.copy()

# Split into train/test using the same approach as the main modeling section
n_total = len(df_sim)
n_test = int(n_total * config.TEST_SIZE)
n_train = n_total - n_test

# Load the feature selection mask for linear models
selected_features_mask = joblib.load("./assets/saved_results/selected_features_mask.pkl")
linear_feature_cols = [col for i, col in enumerate(feature_cols) if selected_features_mask[i]]
tree_feature_cols = feature_cols  # Use all features for tree models

# Prepare full training data (will be expanded in the adaptive strategy)
X_train_full_sim = df_sim.iloc[:n_train]
y_train_full_sim = df_sim.iloc[:n_train]

# Prepare the 180-day test slice for the simulation
X_test_sim = df_sim.iloc[n_train : n_train + 180]
y_test_sim = df_sim.iloc[n_train : n_train + 180]

# Pre-slice feature sets for convenience
X_test_sim_linear = X_test_sim[linear_feature_cols]
X_test_sim_tree = X_test_sim[tree_feature_cols]

print(f"Initial training data size: {len(X_train_full_sim)} days")
print(f"Simulation test data size: {len(X_test_sim)} days")

Initial training data size: 3310 days
Simulation test data size: 180 days

static_results = dd.run_static_prediction(
    ensemble_models=sim_ensemble,
    X_test_linear=X_test_sim_linear,
    X_test_tree=X_test_sim_tree,
    y_test=y_test_sim,
)

# Store the baseline performance metrics
static_predictions = static_results["predictions"]
static_rmse = static_results["avg_metrics"]["rmse"]
static_r2 = static_results["avg_metrics"]["r2"]

Making predictions for all 5 horizons (Static Baseline)...

  Predicting horizon h1 (t+1 days ahead)...
  Predicting horizon h2 (t+2 days ahead)...
  Predicting horizon h3 (t+3 days ahead)...
  Predicting horizon h4 (t+4 days ahead)...
  Predicting horizon h5 (t+5 days ahead)...

============================================================
✓ Static Ensemble Performance by Horizon:
============================================================
Horizon      RMSE (°C)       MAE (°C)        R²        
------------------------------------------------------------
t+1          0.7625          0.6154          0.7152    
t+2          0.9418          0.7491          0.5652    
t+3          1.0482          0.8452          0.4608    
t+4          1.0811          0.8714          0.4264    
t+5          1.0642          0.8612          0.4469    
------------------------------------------------------------
Average      0.9796          0.7885          0.5229    
============================================================

# --- 1. Initialize the Drift Detector ---
detector = dd.EnsembleDriftDetector(baseline_rmse=static_rmse, baseline_r2=static_r2)
# Set the training data distributions as the reference
detector.set_baseline_distributions(
    X_train_full_sim[tree_feature_cols],
    y_train_full_sim["target_temp_t+1"],
)

# --- 2. Run the Adaptive Simulation ---
# adaptive_results = dd.run_adaptive_simulation(
#     initial_models=sim_ensemble,
#     drift_detector=detector,
#     X_train_full=X_train_full_sim,
#     y_train_full=y_train_full_sim,
#     X_test_full=X_test_sim,
#     y_test_full=y_test_sim,
#     linear_feature_cols=linear_feature_cols,
#     tree_feature_cols=tree_feature_cols,
# )
adaptive_results = joblib.load("./assets/saved_results/q7_adaptive_results.pkl")

# --- 3. Display and Interpret Results ---
adaptive_rmse = adaptive_results["avg_metrics"]["rmse"]
adaptive_r2 = adaptive_results["avg_metrics"]["r2"]
retrain_dates = adaptive_results["retrain_dates"]
retrain_priorities = adaptive_results["retrain_priorities"]

✓ Drift Detector initialized (3-TIER PRIORITY SYSTEM)
  Baseline RMSE: 0.9796°C
  Baseline R²: 0.5229

  Priority Thresholds:
  🔴 HIGH (Immediate):  RMSE>1.5× OR R²drop>0.3
  🟡 MEDIUM (Proactive): Data+Label drift OR moderate performance drop
  🟢 LOW (Scheduled):    90 days without retrain
✓ Stored 139 feature distributions

dd.plot_simulation_results(
    static_results=static_results,
    adaptive_results=adaptive_results,
    y_test=y_test_sim,
)

df_hourly_clean = hp.load_and_clean_hourly_data(config.HOURLY_DATA_PATH)

Raw hourly data loaded: (94248, 25)

Dropped 28 duplicate rows.

✓ Data cleaning complete. No remaining NaNs.

df_s1_features = hp.create_daily_enriched_features(df_hourly_clean)

Daily aggregation complete. Shape: (3927, 51)

✓ Feature engineering complete. Final shape: (3927, 154)

# This function handles data splitting, LASSO, model training, and evaluation.
# results_s1_df = hp.run_strategy_1_evaluation(
#     df_features=df_s1_features,
#     test_size=config.TEST_SIZE
# )
results_s1_df = joblib.load("./assets/saved_results/results_s1_df.pkl")

table = Table(title="[bold]Strategy 1 Results: Aggregate-then-Predict[/bold]")
table.add_column("Model", style="cyan", no_wrap=True)
table.add_column("Avg R² ↑", style="green", justify="right")
table.add_column("Avg RMSE ↓", style="yellow", justify="right")
table.add_column("Avg MAE ↓", style="yellow", justify="right")

for _, row in results_s1_df.iterrows():
    table.add_row(
        row["Model"],
        f"{row['Avg_R2']:.4f}",
        f"{row['Avg_RMSE']:.4f}°C",
        f"{row['Avg_MAE']:.4f}°C"
    )
console.print(table)

      Strategy 1 Results: Aggregate-then-Predict      
┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━┓
┃ Model          ┃ Avg R² ↑ ┃ Avg RMSE ↓ ┃ Avg MAE ↓ ┃
┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━┩
│ Champion_Daily │   0.6195 │   0.9664°C │  0.7698°C │
│ Q8.1_Ensemble  │   0.6117 │   0.9762°C │  0.7730°C │
│ Q8.1_RidgeCV   │   0.5969 │   0.9944°C │  0.7870°C │
└────────────────┴──────────┴────────────┴───────────┘

# This function handles the hourly feature engineering, 120-hour target creation, and data splitting.
s2_data = hp.prepare_strategy_2_data(
    df_hourly_clean=df_hourly_clean,
    test_start_timestamp=pd.to_datetime("2024-02-20"),
)

Hourly feature engineering complete. Shape: (94220, 95)

Data prepared for Strategy 2. Train shape: (79993, 215), Test shape: (14035, 215)

# This function runs the entire S2 pipeline: trains 5 models, predicts, aggregates, and evaluates.
# s2_results = hp.run_strategy_2_evaluation(
#     X_train=s2_data["X_train"],
#     y_train=s2_data["y_train"],
#     X_test=s2_data["X_test"],
#     y_test=s2_data["y_test"],
#     day_targets=s2_data["day_targets"],
# )

# results_s2_df = pd.DataFrame([s2_results])
results_s2_df = joblib.load("./assets/saved_results/results_s2_df.pkl")

final_hourly_results_df = pd.concat([results_s1_df, results_s2_df], ignore_index=True)
final_hourly_results_df = final_hourly_results_df.sort_values(by="Avg_R2", ascending=False)

table = Table(title="[bold]Strategy 2 Results: Predict-then-Aggregate[/bold]")
table.add_column("Model", style="cyan", no_wrap=True)
table.add_column("Avg R² ↑", style="green", justify="right")
table.add_column("Avg RMSE ↓", style="yellow", justify="right")
table.add_column("Avg MAE ↓", style="yellow", justify="right")

for _, row in final_hourly_results_df.iterrows():
    table.add_row(
        row["Model"],
        f"{row['Avg_R2']:.4f}",
        f"{row['Avg_RMSE']:.4f}°C",
        f"{row['Avg_MAE']:.4f}°C"
    )
console.print(table)

          Strategy 2 Results: Predict-then-Aggregate          
┏━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━┓
┃ Model                  ┃ Avg R² ↑ ┃ Avg RMSE ↓ ┃ Avg MAE ↓ ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━┩
│ Champion_Daily         │   0.6195 │   0.9664°C │  0.7698°C │
│ Q8.2_Hourly_MultiModel │   0.6185 │   0.9436°C │  0.7569°C │
│ Q8.1_Ensemble          │   0.6117 │   0.9762°C │  0.7730°C │
│ Q8.1_RidgeCV           │   0.5969 │   0.9944°C │  0.7870°C │
└────────────────────────┴──────────┴────────────┴───────────┘

hourly_analysis_results = hp.run_hourly_deep_dive_analysis(
    X_test=s2_data["X_test"],
    y_test=s2_data["y_test"],
    day_targets=s2_data["day_targets"],
    model_path_template="./assets/models/hourly/hourly_model_day_{day_num}.pkl",
)

--- Overall Average Hourly Performance ---
  - Avg. Hourly RMSE: 1.5465°C
  - Avg. Hourly R²:   0.7282

ridge_model = joblib.load("./assets/models/daily/RidgeCV_h1.pkl")
data_input_dir = "./assets/saved_results/step_9_data"
file_X_test_linear = os.path.join(data_input_dir, "X_test_linear.joblib")
file_X_test_tree = os.path.join(data_input_dir, "X_test_tree.joblib")

X_test_linear = joblib.load(file_X_test_linear)
X_test_tree = joblib.load(file_X_test_tree)

print(f"Shape of X_test_linear: {X_test_linear.shape}")
print(f"Shape of X_test_tree: {X_test_tree.shape}")

Shape of X_test_linear: (584, 41)
Shape of X_test_tree: (584, 139)

# --- 1. Define the input format for the model ---
initial_type_ridge = [("float_input", FloatTensorType([None, X_test_linear.shape[1]]))]

# ---2. Convert the Scikit-learn model to the ONNX format ---
onnx_model_ridge = skl2onnx.convert_sklearn(
    ridge_model,
    initial_types=initial_type_ridge,
    target_opset=12,  # Stable compatibility
)

# --- 3. Save the ONNX model to a file ---
onnx_filename_ridge = "./assets/models/onnx/ridge_forecaster.onnx"
with open(onnx_filename_ridge, "wb") as f:
    f.write(onnx_model_ridge.SerializeToString())

benchmark_results = {}

# Prepare test data (ONNX Runtime requires float32)
X_test_linear_np = X_test_linear.astype(np.float32).values

# Load the ONNX Runtime session for the converted RidgeCV model.
sess_ridge = rt.InferenceSession(onnx_filename_ridge)
input_name_ridge = sess_ridge.get_inputs()[0].name
output_name_ridge = sess_ridge.get_outputs()[0].name

# Run prediction with the original Scikit-learn model.
pred_original_ridge = ridge_model.predict(X_test_linear_np)

# Run prediction with the ONNX model.
pred_onnx_ridge = sess_ridge.run([output_name_ridge], {input_name_ridge: X_test_linear_np})[0]

# Compare predictions with numerical tolerance
np.testing.assert_allclose(pred_original_ridge, pred_onnx_ridge.flatten(), rtol=1e-5, atol=1e-5)

# Use %timeit with the '-o' flag to capture the timing result as an object.
print("Timing original Scikit-learn model (.pkl):")
t_original_ridge = %timeit -n 100 -r 7 -o ridge_model.predict(X_test_linear_np)
benchmark_results["RidgeCV (.pkl)"] = t_original_ridge

print("\nTiming ONNX Runtime model (.onnx):")
t_onnx_ridge = %timeit -n 100 -r 7 -o sess_ridge.run([output_name_ridge], {input_name_ridge: X_test_linear_np})
benchmark_results["RidgeCV (ONNX)"] = t_onnx_ridge


processed_data = []
for name, result in benchmark_results.items():
    model_name, model_type = name.split(" ")

    # Calculate throughput (predictions per second)
    throughput = len(X_test_linear_np) / result.average
    processed_data.append({
        "Model": model_name,
        "Deployment Type": model_type.strip("()"),
        # result.average is in seconds; convert to microseconds (µs) for readability.
        "Avg. Time (µs)": result.average * 1_000_000,
        "Std. Dev. (µs)": result.stdev * 1_000_000,
        "Throughput (predictions/sec)": throughput,
    })

df_results = pd.DataFrame(processed_data)

# Calculate the speedup factor relative to the .pkl baseline.
baseline_time = df_results.loc[0, "Avg. Time (µs)"]
onnx_time = df_results.loc[1, "Avg. Time (µs)"]
df_results.loc[1, "Speedup"] = baseline_time / onnx_time
df_results.loc[0, "Speedup"] = 1.0


result_table = (
    df_results.style.format({
        "Avg. Time (µs)": "{:,.1f}".format,
        "Std. Dev. (µs)": "± {:,.1f}".format,
        "Speedup": "{:,.2f}x".format,
        "Throughput (predictions/sec)": "{:,.0f}".format,
    })
    .set_properties(**{"text-align": "center"})
    .set_table_styles([dict(selector="th", props=[("text-align", "center")])])
)

print("--- Deployment Performance ---")
display(result_table)

Timing original Scikit-learn model (.pkl):
249 μs ± 40.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Timing ONNX Runtime model (.onnx):
57.6 μs ± 8.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
--- Deployment Performance ---

"./ui/lasso_sunburst_chart.html"

'./ui/lasso_sunburst_chart.html'

	count	mean	std	min	25%	50%	75%	max
latitude	3927.0	10.776000	0.000000e+00	10.776	10.776	10.776	10.776	10.776
longitude	3927.0	106.701000	1.421266e-14	106.701	106.701	106.701	106.701	106.701
tempmax	3927.0	33.071683	1.791789e+00	24.600	32.000	33.000	34.000	39.000
tempmin	3927.0	25.136949	1.586575e+00	18.000	24.000	25.000	26.000	30.000
temp	3927.0	28.445837	1.386448e+00	22.800	27.500	28.400	29.300	33.000
feelslikemax	3927.0	38.463305	3.247412e+00	24.600	36.400	38.800	40.600	48.700
feelslikemin	3927.0	25.830736	2.943036e+00	18.000	24.000	25.000	26.000	40.000
feelslike	3927.0	31.671760	2.958532e+00	22.800	29.700	31.400	33.500	42.000
dew	3927.0	23.503896	2.242134e+00	12.700	22.400	24.200	25.000	27.900
humidity	3927.0	76.570512	9.681565e+00	49.500	69.700	77.500	84.200	99.400
precip	3927.0	5.208251	1.237959e+01	0.000	0.100	0.800	5.000	227.200
precipprob	3927.0	75.986758	4.272181e+01	0.000	100.000	100.000	100.000	100.000
precipcover	3927.0	10.342480	1.321839e+01	0.000	4.170	4.170	12.500	100.000
windgust	3927.0	31.897428	1.051264e+01	8.300	24.800	29.900	36.700	216.000
windspeed	3927.0	19.462210	5.504470e+00	7.600	15.800	18.400	22.300	50.000
windspeedmax	3927.0	19.462210	5.504470e+00	7.600	15.800	18.400	22.300	50.000
windspeedmean	3927.0	10.029870	3.400972e+00	3.200	7.400	9.400	12.300	23.000
windspeedmin	3927.0	3.253603	2.866895e+00	0.000	1.400	2.300	4.900	15.700
winddir	3927.0	185.397479	8.071268e+01	0.000	127.750	180.900	251.600	359.700
sealevelpressure	3927.0	1009.114770	2.094073e+00	1001.500	1007.700	1009.000	1010.500	1016.200
cloudcover	3927.0	55.594321	1.359861e+01	9.700	46.450	55.000	64.800	100.000
visibility	3927.0	9.614719	7.783575e-01	5.400	9.200	9.800	10.200	20.000
solarradiation	3927.0	213.147670	4.969988e+01	0.000	185.800	221.200	250.000	315.800
solarenergy	3927.0	18.403489	4.296753e+00	0.000	16.100	19.100	21.600	27.200
uvindex	3927.0	7.803667	1.494908e+00	0.000	7.000	8.000	9.000	10.000
severerisk	1196.0	27.604515	1.737145e+01	5.000	10.000	30.000	30.000	75.000
moonphase	3927.0	0.483990	2.887589e-01	0.000	0.250	0.500	0.750	0.980

Feature Group	Feature Name	Data Type	Value Range (numerical) \| Example Values (categorical)	Definition & How it may helps predicting target variable `temp`
Identifier/Metadata	`name`, `address`...	`object`	`['Hồ Chí Minh city']`	Constant metadata identifying the location. Action: Will be dropped.
	`datetime`	`object`	`2015-01-01` to `2025-10-01`	The unique timestamp for each daily record. Action: Will be converted to the DataFrame index.
Target & Temp.	`temp`	`float64`	`22.8` - `33.0` (°C)	Primary Target Variable. The average daily temperature.
	`tempmax`/`tempmin`	`float64`	`24.6`-`39.0` / `18.0`-`30.0` (°C)	The daily maximum and minimum temperatures. Key drivers of the daily average.
	`feelslike`...	`float64`	e.g., `22.8` - `42.0` (°C)	The perceived temperature, accounting for humidity and wind. Highly correlated with `temp`.
Atmospheric	`humidity`	`float64`	`49.5` - `99.4` (%)	Relative humidity. Crucial in a tropical climate as high humidity traps heat, especially impacting overnight lows.
	`dew`	`float64`	`12.7` - `27.9` (°C)	Dew Point. The temperature at which air becomes saturated. A direct measure of atmospheric moisture.
	`sealevelpressure`	`float64`	`1001.5` - `1016.2` (hPa)	Air pressure adjusted to sea level. Changes can indicate shifting weather patterns (e.g., low pressure is associated with storms).
	`cloudcover`	`float64`	`9.7` - `100.0` (%)	The percentage of the sky covered by clouds. Directly impacts solar radiation reaching the surface.
	`visibility`	`float64`	`5.4` - `20.0` (km)	The distance at which objects can be clearly seen. Often reduced by rain or high humidity.
Precipitation	`precip`	`float64`	`0.0` - `227.2` (mm)	Total daily precipitation. Defines the wet vs. dry seasons, which are central to HCMC's climate.
	`precipprob`	`int64`	`0` or `100` (%)	Probability of precipitation. In this dataset, there's only 2 values, can be understood as 0 or 1.
	`preciptype`	`object`	`['rain', nan]`	The type of precipitation. `NaN` values deterministically correspond to `precip = 0`. Action: Impute `NaN` with `'none'`.
Wind	`windspeed`	`float64`	`7.6` - `50.0` (km/h)	The maximum sustained wind speed for the day.
	`windgust`	`float64`	`8.3` - `216.0` (km/h)	The highest instantaneous wind speed recorded. Can indicate storm activity.
	`winddir`	`float64`	`0.0` - `359.7` (°)	The direction from which the wind originates. Seasonal wind patterns (monsoons) are key climate drivers.
Solar & Celestial	`solarradiation`	`float64`	`0.0` - `315.8` (W/m²)	The amount of solar energy reaching the surface. The primary driver of daytime temperature.
	`solarenergy`	`float64`	`0.0` - `27.2` (MJ/m²)	The total solar energy over the day. Highly correlated with `solarradiation`.
	`uvindex`	`int64`	`0` - `10`	The strength of ultraviolet radiation. A proxy for sun intensity and lack of cloud cover.
	`sunrise`/`sunset`	`object`	`2015-01-01T06:11:22`	Timestamps for sunrise and sunset. Action: Can be engineered into a `daylight_hours` feature.
	`moonphase`	`float64`	`0.0` (new) - `0.98` (waning)	The fraction of the moon that is illuminated. May have a subtle, long-term cyclical influence.
Text & Categorical	`conditions`	`object`	`['Partially cloudy', 'Rain...']`	A high-level summary of the day's weather.
	`description`	`object`	`['Partly cloudy throughout...']`	A more detailed, human-readable description of the weather. Action: Will be processed using NLP techniques.
	`icon`	`object`	`['partly-cloudy-day', 'rain']`	A machine-readable categorical summary of conditions.
Risk	`severerisk`	`float64`	`5.0` - `75.0`	A numerical rating for the risk of severe weather events. `NaN` indicates no recorded risk. Action: Impute `NaN` with `0`.

	year	month	week_of_year	day_of_year	sin_day_of_year	cos_day_of_year	sin_month	cos_month	fourier_sin_2	fourier_cos_2	fourier_sin_3	fourier_cos_3
datetime
2015-01-01	2015	1	1	1	0.017166	0.999853	0.5	0.866025	0.034328	0.999411	0.051479	0.998674
2015-01-02	2015	1	1	2	0.034328	0.999411	0.5	0.866025	0.068615	0.997643	0.102821	0.994700
2015-01-03	2015	1	1	3	0.051479	0.998674	0.5	0.866025	0.102821	0.994700	0.153891	0.988088
2015-01-04	2015	1	1	4	0.068615	0.997643	0.5	0.866025	0.136906	0.990584	0.204552	0.978856
2015-01-05	2015	1	2	5	0.085731	0.996318	0.5	0.866025	0.170830	0.985301	0.254671	0.967028

	daylight_hours	temp_range	dew_point_deficit
datetime
2015-01-01	11.504722	8.0	8.7
2015-01-02	11.507222	10.0	9.2
2015-01-03	11.509722	9.0	7.8
2015-01-04	11.512500	8.0	5.8
2015-01-05	11.515556	5.7	4.6

	text_svd_0	text_svd_1	text_svd_2	text_svd_3	text_svd_4	text_svd_5	text_svd_6
datetime
2015-01-01	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-02	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-03	0.555798	-0.049113	0.585367	0.280771	-0.240929	0.061356	0.448770
2015-01-04	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-05	0.407113	-0.081170	0.087575	0.089282	-0.493769	-0.186068	-0.152486

	tempmax	tempmin	temp	feelslikemax	feelslikemin	feelslike	dew	humidity	precip	precipprob	...	precip_roll_28d_sum	preciptype_none	preciptype_rain	text_svd_0	text_svd_1	text_svd_2	text_svd_3	text_svd_4	text_svd_5	text_svd_6
datetime
2015-01-01	31.0	23.0	26.6	31.4	23.0	26.8	17.9	60.3	0.0	0	...	NaN	1	0	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-02	30.0	20.0	25.0	30.4	20.0	25.1	15.8	57.2	0.0	0	...	0.0	1	0	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-03	32.0	23.0	26.8	33.5	23.0	27.4	19.0	63.7	0.2	100	...	0.0	0	1	0.555798	-0.049113	0.585367	0.280771	-0.240929	0.061356	0.448770
2015-01-04	32.0	24.0	27.1	34.8	24.0	28.3	21.3	71.7	0.0	0	...	0.2	1	0	0.879053	-0.176581	-0.320034	-0.237530	-0.112032	-0.016034	0.155117
2015-01-05	30.6	24.9	26.7	33.2	24.9	27.7	22.1	76.4	6.0	100	...	0.2	0	1	0.407113	-0.081170	0.087575	0.089282	-0.493769	-0.186068	-0.152486

METRIC	FORMULA	CORE QUESTION	INTEPRETATIONS
R-squared ($R^2$)	$1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2}$	"How much of the temperature's variance is explained by our model?"	A scale-free score from $(-\infty, 1]$. A score of 0.7 means our model explains 70% of the variance. Primary metric for model comparison.
Adjusted $R^2$	$1 - \frac{(1-R^2)(n-1)}{n-p-1}$	"What is the $R^2$, after penalizing for adding more features?"	Useful for explanatory linear modeling to avoid overfitting by adding useless predictors. Not used in our final reports.
RMSE	$\sqrt{\frac{1}{n} \sum_{i=1}^{n}(y_i - \hat{y}_i)^2}$	"What is the typical magnitude of our error, penalizing larger errors more?"	Also in degrees Celsius. Because errors are squared, RMSE is more sensitive to large forecast misses than MAE. A key indicator of error size.
MAE	$\frac{1}{n} \sum_{i=1}^{n} \\|y_i - \hat{y}_i\\|$	"On average, what is the absolute error of our forecast in degrees?"	Highly interpretable. An MAE of 0.5 means our forecasts are off by an average of 0.5°C. Robust to outliers.
MAPE	$\frac{100}{n} \sum_{i=1}^{n} \\|\frac{y_i - \hat{y}_i}{y_i}\\|$	"On average, what is the percentage error of our forecast?"	A relative error metric. A MAPE of 2% means our forecasts are off by an average of 2% of the actual temperature.

Priority	Trigger Condition	Thresholds (from our simulation)	Action
🔴 HIGH	Severe Performance Degradation: A sudden, unacceptable drop in accuracy.	`RMSE Ratio > 1.50` OR `R² Drop > 0.30`	Immediate Retrain
🟡 MEDIUM	Moderate Performance Drop OR Combined Distributional Drift	`RMSE Ratio > 1.30` OR `R² Drop > 0.20` OR (`Data Drift` AND `Label Drift`)	Proactive Retrain
🟢 LOW	Scheduled Maintenance: The model has not been retrained for an extended period.	`Days since last retrain ≥ 90`	Scheduled Retrain

Challenge	ONNX Advantage
Interoperability	Convert both RidgeCV (Scikit-learn) and LGBMRegressor (LightGBM) into `.onnx` format → deployment only requires ONNX Runtime, not multiple library dependencies.
Performance	ONNX Runtime, written in C++, supports hardware acceleration (e.g., AVX, GPU), enabling significantly faster inference compared to using `.predict()` in Python.
Cross-Platform Deployment	Although training models in Python, we can deploy them in C#, C++, Java, or JavaScript environments. Easily deploy models to edge devices, mobile, and IoT systems with hardware constraints, or to cloud platforms (Azure ML, AWS, GCP) using the standardized ONNX format.
Production Performance & Scalability	Ideal for real-time inference where high speed and low latency are required; ONNX Runtime optimizes for production-grade workloads.
Team Collaboration & Framework Independence	Different teams using different frameworks can share models through the ONNX standard. This avoids vendor lock-in and allows flexible switching between frameworks and tools.
Model Optimization	Tools like quantization reduce model size and improve inference speed — ideal for edge and resource-limited devices.

	Model	Deployment Type	Avg. Time (µs)	Std. Dev. (µs)	Throughput (predictions/sec)	Speedup
0	RidgeCV	.pkl	249.4	± 40.3	2,342,028	1.00x
1	RidgeCV	ONNX	57.6	± 8.1	10,144,122	4.33x

Ho Chi Minh City Temperature Forecasting Group 2 - Analysis & Implementation

0. Setup & Global Configuration

0. Setup & Global Configuration

0.1. Centralized Pipeline Configuration

0.1. Centralized Pipeline Configuration

Q1. Data Sourcing

Q1. Data Sourcing

Q2. Data Understanding & Exploratory Data Analysis (EDA)

Q2. Data Understanding & Exploratory Data Analysis (EDA)

Q2.1. Initial Data Load & Quality Audit

Q2.1. Initial Data Load & Quality Audit

Q2.2. Feature Definitions & Domain Context

Q2.2. Feature Definitions & Domain Context

Q2.3. Deep Dive: Key Predictive Feature Groups

Q2.3. Deep Dive: Key Predictive Feature Groups

The Solar Radiation & Cloud Cover Nexus¶

The Atmospheric Moisture Duo: Humidity & Dew Point¶

The Precipitation System: A Cooling & Lagged Effect¶

Wind Dynamics: Direction & Speed¶

Q2.4. Target Variable Analysis: Trends & Seasonality

Q2.4. Target Variable Analysis: Trends & Seasonality

Q2.5. Seasonal Pattern Deep Dive: The Monsoon Cycle

Q2.5. Seasonal Pattern Deep Dive: The Monsoon Cycle

Q2.6. STL Decomposition: Isolating Key Signals

Q2.6. STL Decomposition: Isolating Key Signals

Q2.7. Predictor Analysis: Distribution, Skew, & Outliers

Q2.7. Predictor Analysis: Distribution, Skew, & Outliers

Q2.8. Bivariate & Seasonal Relationship Analysis

Q2.8. Bivariate & Seasonal Relationship Analysis

Q3. Data Processing & Correlation Analysis

Q3. Data Processing & Correlation Analysis

Q3.1. Column Sanitization

Q3.1. Column Sanitization

Q3.2. Missing Value Imputation

Q3.2. Missing Value Imputation

Q3.3. Correlations Check

Q3.3. Correlations Check

Q3.4. Strategy for Multicollinearity

Q3.4. Strategy for Multicollinearity

Q3.5. Final Processed DataFrame before Feature Engineering

Q3.5. Final Processed DataFrame before Feature Engineering

Q4. Feature Engineering

Q4. Feature Engineering

Q4.1. Foundational Temporal Features

Q4.1. Foundational Temporal Features

Q4.2. Domain-Specific Feature Creation

Q4.2. Domain-Specific Feature Creation

Q4.3. Lag & Rolling Window Features

Q4.3. Lag & Rolling Window Features

Q4.4. Categorical Feature Encoding

Q4.4. Categorical Feature Encoding

Q4.5. Unstructured Text Feature Engineering

Q4.5. Unstructured Text Feature Engineering

Q4.6. Final Verification & NaN Handling

Q4.6. Final Verification & NaN Handling

Q5. Modeling Strategy & Evaluation

Q5. Modeling Strategy & Evaluation

Q5.1. Forecasting Strategies, Target Definition & Data Splitting

Q5.1. Forecasting Strategies, Target Definition & Data Splitting

Q5.1.1. Defining the Forecasting Strategy

Q5.1.1. Defining the Forecasting Strategy

Forecasting Architectures¶

What to Predict?¶

Q5.1.2. Target Variable Creation & Final Alignment

Q5.1.2. Target Variable Creation & Final Alignment

Q5.1.3. The Final Train-Test Split

Q5.1.3. The Final Train-Test Split

Q5.1.4. Creating Specialized Feature Subsets

Q5.1.4. Creating Specialized Feature Subsets

Q5.2. Evaluation Framework

Q5.2. Evaluation Framework

Q5.2.1. Evaluation Metrics

Q5.2.1. Evaluation Metrics

Q5.2.2. Cross-Validation Strategy

Q5.2.2. Cross-Validation Strategy

Q5.3. The Baseline Model: Establishing the Performance Floor

Q5.3. The Baseline Model: Establishing the Performance Floor

Q5.4. Model Bake-Off

Q5.4. Model Bake-Off

Q5.4.1. Stage 1: High-Level Leaderboard

Ho Chi Minh City Temperature Forecasting
Group 2 - Analysis & Implementation